* This file creates figures verifying robutness of reults to covariate imbalance in the input experiment
* Paper Figures S4-S8

* Variables for entropy balancing
local bal_resp m_gender m_age m_educ_p m_educ_s f_gender f_age f_educ_p f_educ_s
local bal_hh hh_size caste_scst bl_pulse_prev bl_wealth_win5 bl_land_win5

* Make graph pretty
local grset graphregion(color(white)) ylabel(, angle(0)) grid(glcolor(gs0) glpattern(dash) glwidth(thin) nogextend) plotregion(margin(0) lcolor(black)) 

tempfile treats
graph drop _all
snapshot erase _all
local num = 4

*****
* Household treatment assignemnts
*****

use "$out_data/survey_roster.dta", clear
merge m:1 village using "$admin_data/treatment_assignment.dta", keepusing(treat_input treat_ext) assert(match using)
drop if _m==2
drop _m

* Split land and wealth into above/below median for heterogeneity analysis
foreach v in land wealth {
	quietly sum bl_`v'_win5, d
	gen byte bl_`v'_high = (bl_`v'_win5 > `r(p50)') if !missing(bl_`v'_win5)
}

* Gen stuff for entropy balancing
recode caste (1 4 = 0) (2 3 = 1), gen(caste_scst)
recode m_education (1 = 0) (2/7 = 1), gen(m_educ_primary)
recode f_education (1 = 0) (2/7 = 1), gen(f_educ_primary)
recode m_education (1/2 = 0) (3/7 = 1), gen(m_educ_secondary)
recode f_education (1/2 = 0) (3/7 = 1), gen(f_educ_secondary)

ebalfit `bal_resp' `bal_hh', by(treat_input) gen(ebals)

merge m:1 hhid using "$out_data/survey_drops.dta", assert(match) nogen

// Make year*treat interactions
expand 3
bys hhid: gen byte year = _n
* Dummy out treatment interactions
gen byte B_y2 = year==2
label var B_y2 "Year 2"
gen byte B_y3 = year==3
label var B_y3 "Year 3"
gen byte B_y1_t = (year==1)*(treat_input==1)
label var B_y1_t "Yr. 1"
gen byte B_y2_t = (year==2)*(treat_input==1)
label var B_y2_t "Yr. 2"
gen byte B_y3_t = (year==3)*(treat_input==1)
label var B_y3_t "Yr. 3"

gen byte B_y1_s = (year==1)*(treat_ext==1)
label var B_y1_s "Subs. Yr. 1"
gen byte B_y2_s = (year==2)*(treat_ext==1)
label var B_y2_s "Subs. Yr. 2"
gen byte B_y3_s = (year==3)*(treat_ext==1)
label var B_y3_s "Subs. Yr. 3"
gen byte B_y1_e = (year==1)*(treat_ext==2)
label var B_y1_e "Ext. Yr. 1"
gen byte B_y2_e = (year==2)*(treat_ext==2)
label var B_y2_e "Ext. Yr. 2"
gen byte B_y3_e = (year==3)*(treat_ext==2)
label var B_y3_e "Ext. Yr. 3"

* Generate round-specific baseline interactions for double lasso
foreach vbl in $controls {
	if ("`vbl'"=="i.block") continue
	
	* Parse out factor variables
	if substr("`vbl'",2,1)=="." {
		local vname = substr("`vbl'",3,.)
		quietly tab `vname', gen(`vname'_)
		local items = "`vname'_*"
	}
	else local items = "`vbl'"
	
	foreach item of varlist `items' {
	forvalues yr = 1/3 {
		gen dl_y`yr'_`item' = `item'*(year==`yr')
	}
	}
}

save `treats'

*****
* Adoption and Area; Production; Yield
*****

use "$out_data/outcome_adoption_production.dta", clear
collapse (sum) pulse_prod_win5 pulse_area, by(hhid year season) fast

* Make adoption dummy
gen byte pulse_adopt = pulse_area>0 if !missing(pulse_area)

* Make pulse yield (total)
gen pulse_yield = pulse_prod_win5 / pulse_area

* Reshape over seasons and merge in pulse stocks
reshape wide pulse_adopt pulse_area pulse_prod_win5 pulse_yield, i(hhid year) j(season)
merge 1:1 hhid year using "$out_data/outcome_stocks.dta", assert(match) nogen keepusing(pulse_time*)
egen pulse_time_tot = rowmax(pulse_time_*)

* Offset pulse stocks by a year
sort hhid year
by hhid (year): gen pulse_time_prev = pulse_time_tot[_n+1]

* Add in treatment data
merge 1:1 hhid year using `treats', assert(match) nogen


*****
* Regressions: Adoption
*****
local name adopt
eststo clear

local v_`name' pulse_adopt1 pulse_area1 pulse_adopt2 pulse_area2 pulse_adopt3 pulse_area3
	label var pulse_adopt1 "Kharif Adoption"
	label var pulse_area1 "Kharif Area"
	label var pulse_adopt2 "Rabi Adoption"
	label var pulse_area2 "Rabi Area"
	label var pulse_adopt3 "Zaid Adoption"
	label var pulse_area3 "Zaid Area"
snapshot save, label("s_`name'")


foreach vbl of varlist `v_`name'' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy

	* Regressions for figure
	eststo `vbl'_m: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	eststo `vbl'_d: dsregress `vbl' $B_main i.block if drop_`samp'==0, controls(dl_*) vce(cluster village)
	eststo `vbl'_e: reg `vbl' $B_main $controls if drop_`samp'==0 [aweight=ebals], vce(cluster village)
	eststo `vbl'_n: reg `vbl' $B_main i.block if drop_`samp'==0, vce(cluster village)
	
	coefplot `vbl'_m `vbl'_d `vbl'_e `vbl'_n ///
	, keep($B_coef ) order($B_coef ) xline(0) `grset' name(`vbl') title("`: variable label `vbl''") ///
	p1(label("Full Controls") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin))) ///
	p2(label("Double Lasso") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p3(label("Entropy Balance") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p4(label("No Controls") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin)))

	graph export "$out_figures/Parts/fig_controls_`vbl'.png", replace
}

grc1leg `v_`name'', cols(2) graphregion(color(white)) ysize(8)
graph export "$out_figures/Figure_S`num++'.png", replace
graph drop _all

*****
* Regressions: Yield
*****
local name yield
eststo clear

local vbls pulse_yield1 pulse_yield2 pulse_yield3
/// NO PCA: only 43 obs have yield in each season
// pca `vbls'
// predict pca_yield, score
// local v_yield "`vbls' pca_yield"
local v_`name' "`vbls'"
	label var pulse_yield1 "Kharif Yield"
	label var pulse_yield2 "Rabi Yield"
	label var pulse_yield3 "Zaid Yield"
snapshot save, label("s_`name'")

foreach vbl of varlist `v_`name'' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	* Regressions for figure
	eststo `vbl'_m: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	eststo `vbl'_d: dsregress `vbl' $B_main i.block if drop_`samp'==0, controls(dl_*) vce(cluster village)
	eststo `vbl'_e: reg `vbl' $B_main $controls if drop_`samp'==0 [aweight=ebals], vce(cluster village)
	eststo `vbl'_n: reg `vbl' $B_main i.block if drop_`samp'==0, vce(cluster village)
	
	coefplot `vbl'_m `vbl'_d `vbl'_e `vbl'_n ///
	, keep($B_coef ) order($B_coef ) xline(0) `grset' name(`vbl') title("`: variable label `vbl''") ///
	p1(label("Full Controls") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin))) ///
	p2(label("Double Lasso") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p3(label("Entropy Balance") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p4(label("No Controls") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin)))

	graph export "$out_figures/Parts/fig_controls_`vbl'.png", replace
}

grc1leg `v_`name'', cols(2) graphregion(color(white)) ysize(8)
graph export "$out_figures/Figure_S`num++'.png", replace
graph drop _all

*****
* Regressions: Production and Stocks
*****
local name prod
eststo clear

local vbls pulse_prod_win51 pulse_prod_win52 pulse_prod_win53 
pca `vbls', comp(1)
predict pca_prod, score
local v_`name' "`vbls' pca_prod pulse_time_prev"
	label var pulse_prod_win51 "Kharif Production"
	label var pulse_prod_win52 "Rabi Production"
	label var pulse_prod_win53 "Zaid Production"
	label var pca_prod "Production Index"
	label var pulse_time_prev "Harvest Stock"
snapshot save, label("s_`name'")

foreach vbl of varlist `v_`name'' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	if ("`vbl'"=="pulse_time_prev") local regs "B_y1_t B_y2_t B_y2"
	else local regs "$B_main"

	* Regressions for figure
	eststo `vbl'_m: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	eststo `vbl'_d: dsregress `vbl' `regs' i.block if drop_`samp'==0, controls(dl_*) vce(cluster village)
	eststo `vbl'_e: reg `vbl' $B_main $controls if drop_`samp'==0 [aweight=ebals], vce(cluster village)
	eststo `vbl'_n: reg `vbl' $B_main i.block if drop_`samp'==0, vce(cluster village)
	
	coefplot `vbl'_m `vbl'_d `vbl'_e `vbl'_n ///
	, keep($B_coef ) order($B_coef ) xline(0) `grset' name(`vbl') title("`: variable label `vbl''") ///
	p1(label("Full Controls") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin))) ///
	p2(label("Double Lasso") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p3(label("Entropy Balance") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p4(label("No Controls") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin)))

	graph export "$out_figures/Parts/fig_controls_`vbl'.png", replace
}

grc1leg `v_`name'', cols(2) graphregion(color(white)) ysize(8)
graph export "$out_figures/Figure_S`num++'.png", replace
graph drop _all

*****
* Regressions: Profits
*****
local name earn
* Regressions with combined treatment
* Profits, Production Revenue, Sales Revenue, Cost, Land Owned
use `treats', clear
merge 1:m hhid year using "$out_data/outcome_profits.dta", assert(match) nogen

*** Main regression
eststo clear

local vbls earn_revenue_win5 earn_sales_win5 earn_cost_win5 farm_area_acre_win5
pca `vbls', comp(1)
predict pca_earn, score
local v_`name' "earn_profit_win5 `vbls' pca_earn"
	label var earn_profit_win5 "Farm Profit"
	label var earn_revenue_win5 "Farm Revenue"
	label var earn_sales_win5 "Crop Sales"
	label var earn_cost_win5 "Farm Cost"
	label var farm_area_acre_win5 "Farm Area"
	label var pca_earn "Profit Index"
snapshot save, label("s_`name'")

foreach vbl of varlist `v_`name'' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy

	* Regressions for figure
	eststo `vbl'_m: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	eststo `vbl'_d: dsregress `vbl' $B_main i.block if drop_`samp'==0, controls(dl_*) vce(cluster village)
	eststo `vbl'_e: reg `vbl' $B_main $controls if drop_`samp'==0 [aweight=ebals], vce(cluster village)
	eststo `vbl'_n: reg `vbl' $B_main i.block if drop_`samp'==0, vce(cluster village)
	
	coefplot `vbl'_m `vbl'_d `vbl'_e `vbl'_n ///
	, keep($B_coef ) order($B_coef ) xline(0) `grset' name(`vbl') title("`: variable label `vbl''") ///
	p1(label("Full Controls") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin))) ///
	p2(label("Double Lasso") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p3(label("Entropy Balance") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p4(label("No Controls") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin)))

	graph export "$out_figures/Parts/fig_controls_`vbl'.png", replace
}

grc1leg `v_`name'', cols(2) graphregion(color(white)) ysize(8)
graph export "$out_figures/Figure_S`num++'.png", replace
graph drop _all


*****
* Consumption
***** 
local name consume
* Regressions with combined treatment
* Pulse consumption, Protein consumption, Pulse stocks
use `treats', clear
merge 1:m hhid year using "$out_data/outcome_consumption.dta", assert(match master) nogen
merge m:1 hhid year using "$out_data/outcome_stocks.dta", assert(match) nogen keepusing(stock*)

foreach vbl in kgs kw5 {
	egen stock_`vbl'_tot = rowtotal(stock_`vbl'_*)
}
egen stock_time_tot = rowmax(stock_time_*)

*** Main regression
eststo clear

local vbls stock_kw5_tot stock_time_tot eat_week_pulse_win5_pc eat_day_prot_win5_pc
pca `vbls', comp(1)
pca `vbls', comp(2)
predict pca_eat pca_alt_eat, score
local v_`name' "`vbls' eat_day_pfem_win5 pca_eat"
	label var stock_kw5_tot "Pulse Stock (Kgs.)"
	label var stock_time_tot "Pulse Stock (Mos.)"
	label var eat_week_pulse_win5_pc "Consumption"
	label var eat_day_prot_win5_pc "Daily HH Protein"
	label var eat_day_pfem_win5 "Female Protein"
	label var pca_eat "Consumption Index"
snapshot save, label("s_`name'")

foreach vbl of varlist `v_`name'' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	if ("`vbl'"=="eat_day_pfem_win5") local regs "B_y2_t B_y3_t B_y3"
	else local regs "$B_main"
	
	* Regressions for figure
	eststo `vbl'_m: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	eststo `vbl'_d: dsregress `vbl' `regs' i.block if drop_`samp'==0, controls(dl_*) vce(cluster village)
	eststo `vbl'_e: reg `vbl' $B_main $controls if drop_`samp'==0 [aweight=ebals], vce(cluster village)
	eststo `vbl'_n: reg `vbl' $B_main i.block if drop_`samp'==0, vce(cluster village)
	
	coefplot `vbl'_m `vbl'_d `vbl'_e `vbl'_n ///
	, keep($B_coef ) order($B_coef ) xline(0) `grset' name(`vbl') title("`: variable label `vbl''") ///
	p1(label("Full Controls") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin))) ///
	p2(label("Double Lasso") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p3(label("Entropy Balance") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p4(label("No Controls") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin)))
	
	graph export "$out_figures/Parts/fig_controls_`vbl'.png", replace
}

grc1leg `v_`name'', cols(2) graphregion(color(white)) ysize(8)
graph export "$out_figures/Figure_S`num++'.png", replace
graph drop _all